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I.  INTRODUCTION 


A.  BACKGROUND 

Space  based  operations  pose  significant  operational  and  technical  challenges.  The 
job  of  maintaining  orbit  in  the  space  environment  is  complex,  yet  operators  must  also 
contend  with  outside  threats  to  their  spacecraft.  Threats  to  spacecraft  include  the  space 
medium  itself,  conventional  weapons,  directed  energy  weapons,  and  electronic  warfare 
(Baines,  2006,  pp.  33-39).  One  threat  in  particular  is  growing  at  an  alarming  rate;  the 
threat  of  colliding  with  other  orbiting  objects.  The  advent  of  space  use  has  introduced 
more  than  39,000  traceable  manmade  objects  into  orbit,  with  over  16,000  large  enough  to 
be  currently  tracked  (United  States  Strategic  Command,  2012).  Of  these  objects, 
approximately  five  percent  are  functional  (United  States  Strategic  Command,  2012). 
This  leaves  a  vast  majority  of  objects  orbiting  earth  that  have  no  means  of  maneuvering 
away  from  a  collision. 

This  population  of  orbiting  objects,  including  spacecraft,  rocket  bodies,  and  debris 
is  continuously  growing.  The  risk  of  collision  in  space  rises  with  the  number  of  objects 
in  orbit,  and  has  the  potential  to  create  cascading,  exponential  increases  in  the  number  of 
objects  orbiting  earth  (Kessler  &  Cour-Palais,  1978).  This  increase  is  best  illustrated  by 
the  first  ever  accidental  collision  between  two  satellites,  Iridium  33  and  Cosmos  2251, 
which  left  distinct  shells  of  debris  across  wide  orbits  (Satellite  collision  leaves  significant 
debris  cloud,  2009).  This  collision  alone  produced  more  than  2,000  traceable  objects 
(International  Space  Station  again  dodges  debris,  2011). 

The  situation  is  made  worse  by  the  actions  of  the  international  community.  In 
2007,  China  successfully  tested  a  kinetic  kill  vehicle  against  a  defunct  weather  satellite, 
Fengyun-lC,  creating  approximately  950  traceable  objects  and  many  smaller  objects 
(U.S.  Library  of  Congress,  2007).  This  debris  is  generally  considered  responsible  for  the 
sudden  breakup  of  a  Russian  retro  flection  satellite  in  2013,  which  was  destroyed  after 
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passing  through  the  debris  cloud  left  by  the  Chinese  test  (Kelso,  2013).  As  the 
opportunity  for  collision  increases  due  to  advancing  technology  and  increasing  debris, 
protecting  operational  spacecraft  must  become  a  priority. 

B.  MOTIVATION 

Expanding  reliance  upon  space  assets  and  space  based  capabilities  increases  the 
opportunity  for  and  consequences  of  a  disruption.  Weapon  technology  will  continue  to 
advance  and  spread.  Orbital  debris  can  take  many  years  to  deorbit  and  only  seconds  to 
create.  These  things  will  not  change.  Clearly  satellites  require  greater  levels  of 
protection  than  ever  before,  and  a  promising  technology  may  hold  some  answers  for 
flight  within  and  outside  of  the  atmosphere. 

In  2008,  the  Rockwell  Collins  Company,  sponsored  by  the  Defense  Advanced 
Research  Projects  Agency  (DARPA)  demonstrated  damage  tolerant  controls  (DTC)  in  an 
aircraft  scale  model,  by  jettisoning  60  percent  of  one  wing  and  landing  the  model,  all 
autonomously  (Rockwell  Collins,  2011a).  This  research,  driven  largely  by  increased  the 
United  States  reliance  upon  remotely  piloted  and  autonomous  aircraft,  promises  to 
increase  survivability  and  counteract  the  input  latency  inherent  to  operating  an  unmanned 
aerial  vehicle  (UAV).  DTC  recovers  operational  control  of  a  UAV  before  the  operator 
knows  it  is  damaged,  making  DTC  a  tremendous  asset.  Then,  in  2010,  Rockwell  Collins 
and  DARPA  again  demonstrated  the  ability  to  regain  control  of  a  damaged  autonomous 
UAV,  this  time  the  operational  RQ-7B  Shadow,  continue  the  mission,  and  land 
successfully  (Rockwell  Collins,  2011b).  All  capabilities  were  performed  via  the 
remaining  control  surfaces  only.  This  was  accomplished  using  only  autonomous 
software,  with  modifications  called  adaptive  controls. 

Adaptive  control  is  best  defined  as  “...an  approach  to  dealing  with  uncertain 
systems  or  time-varying  systems”  (Slotine  &  Li,  1991,  p.  204).  This  infers  the  use  of 
what  is  known  as  an  adaptive  controller,  defined  as  “...a  controller  with  adjustable 
parameters  and  a  mechanism  for  adjusting  the  parameters”  (Astrom  &  Wittenmark,  2008, 
p.  1).  By  these  definitions,  adaptive  control  is  a  specific  kind  of  learning  system  that  is 
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generally  considered  appropriate  when  a  system  exhibits  time  variable  dynamics 
(Dumont  &  Huzmean,  2002).  The  demonstration  of  DTC  by  the  RQ-7B  is  enabled  via 
the  use  of  adaptive  control  through  an  adaptive  controller. 

Concurrently,  with  the  successes  of  Rockwell  Collins  and  DARPA,  researchers  at 
the  University  of  Illinois  successfully  modified  the  programming  for  a  UAV  using 
adaptive  controls  (Hovakimyan,  2011).  Using  this  modification  the  UAV  was  able  to 
search  for,  and  follow  a  moving  ground  based  target  autonomously.  These 
demonstrations  represent  industry  firsts  and  a  tremendous  step  forward  in  the  use  of 
adaptive  controls. 

The  successful  tests  by  the  University  of  Illinois,  Rockwell  Collins,  and  DARPA 
with  autonomous  flight  and  adaptive  controls  have  certain  industry  parallels.  Machine 
learning  has  developed  many  optimized  adaptive  controls,  such  as  the  relatively  new 
reward-weighted  regression  model  (Peters  &  Schaal,  2007).  This  model  was  developed 
specifically  to  minimize  the  processing  power  required,  possibly  allowing  for  greater 
adoption  of  autonomous  machine  learning.  Adaptive  controls  have  also  been  designed  to 
adjust  for  a  rapidly  shifting  center  of  gravity  on  light  cargo  trucks  by  the  Robert  Bosch 
company,  dramatically  increasing  stability  and  reducing  rollover  risk  (Liebemann,  Fuher, 
&  Kroger,  2007).  Additionally,  adaptive  controls  are  making  their  way  into  consumer 
products  in  the  fonn  of  vehicular  adaptive  cruise  control  (Man  Truck  &  Bus,  2012). 
Adaptive  cruise  control  is  performed  via  the  controller  maintaining  the  user  set  speed  of 
travel  and  changing  to  match  the  speed  of  slower  moving  vehicles  when  encountered 
(Pananurak,  Somphong,  &  Manukid,  2008).  All  of  these  examples  utilize  adaptive 
controls  to  compensate  for  rapidly  changing  systems. 

Industry  success  and  product  availability  lends  plausibility  to  the  attempt  to  utilize 

adaptive  controls  for  space  based  operations.  The  advantages  of  adaptive  controls  for 

spacecraft  are  seemingly  boundless.  In  a  study  by  the  University  of  Florida,  adaptive 

controls  are  shown  to  overcome  the  real  world  variations  in  torque  seen  when  spacecraft 

utilize  control  movement  gyroscope  (CMG)  gimbals  (MacKunis,  Dixon,  Dupree,  &  Fitz- 

Coy,  2008).  These  variations  in  torque  can  make  spacecraft  attitude  control  all  but 

impossible,  especially  for  small  satellites.  Another  application  is  the  proposed  space  tug, 
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which  moves  between  defunct  orbiting  objects,  collects  them,  and  deorbit  them  (Tewari, 
2012).  Adaptive  controls  would  enable  the  tug  to  compensate  for  unknown  debris  mass 
and  maximize  its  fuel  capabilities. 

As  the  space  medium  becomes  more  crowded  across  all  orbits,  the  development 
of  DTC  for  spacecraft  will  become  a  necessity.  Adaptive  autonomous  operations  has 
been  called  for  in  order  to  improve  spacecraft  survivability  (Baines,  2006,  pp.  45-47). 
Threats  to  spacecraft  exist  at  every  period  of  their  lifetime;  from  launch  to  operations  and 
disposal.  DTC  and  the  recent  developmental  leaps  in  implementing  adaptive  controls 
present  a  possible  answer  to  these  evolving  threats.  This  thesis  will  discuss  the  history, 
problems  associated  with,  and  a  simulated  solution  for  providing  damage  tolerance  via 
adaptive  controls. 

C.  PURPOSE  AND  RESEARCH  QUESTION 

The  purpose  of  this  research  is  to  simulate  a  satellite  with  sudden  mass  loss  and 
simulate  adaptive  controls  to  counteract  the  loss.  DTC  would  be  able  to  react  to  the  loss 
much  quicker  than  any  human  operator,  giving  the  satellite  the  best  opportunity  to 
recover  from  damage.  The  question  this  research  seeks  to  answer  is:  Are  DTC  possible 
for  satellites  in  orbit? 

D.  SCOPE  OF  METHODOLOGY 

The  methodology  of  this  research  begins  with  a  literature  review  of  the 
development  of  adaptive  controls.  Following  that  is  a  detailed  description  of  the 
development  of  the  satellite  simulation.  Last  is  the  data  collection  and  analysis  of  results. 
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II.  LITERATURE  REVIEW 


A.  HISTORICAL  DEVELOPMENT 

In  the  1950s  and  1960s,  the  National  Aeronautics  and  Space  Administration 
(NASA),  the  United  States  Air  Force  (USAF),  and  the  North  American  X-15  rocket 
powered  aircraft  were  pioneering  the  concepts  and  principles  that  would  come  to  define 
modern  powered  flight.  Launched  in  flight  from  the  wing  of  a  B-52A  aircraft,  the  X-15 
would  fire  its  rocket  engine,  achieve  the  mission’s  required  altitude,  and  glide  to  earth  for 
landing.  Among  the  ground  breaking  ideas  proposed  was  a  system  of  adaptive  controls, 
or  a  controller  that  would  take  into  consideration  the  changing  operational  environment  to 
deliver  appropriate  control  to  the  operator  (Staff  of  the  Flight  Research  Center,  1961). 
Primarily  what  the  X-15  proved  was  the  utility  of  something  called  Gain  Scheduling. 

Gain  Scheduling  is  a  method  of  adapting  known  linear  control  technique  to  meet 
the  challenges  required  of  a  nonlinear  system  (Slotine  &  Li,  1991).  By  selecting  enough 
points  across  the  changing  system  and  designing  linear  controls  for  each  portion  of  the 
operation,  it  becomes  possible  to  achieve  adaptive  nonlinear  control.  The  primary 
concern  when  using  gain  scheduling  is  the  lack  of  stability  that  can  occur  if  the  system 
becomes  over  saturated  or  departs  from  the  programmed  model.  In  1967,  oversaturated 
controls  of  the  X-15  that  Major  Mike  Adams  was  piloting  prevented  him  from  regaining 
control  of  the  aircraft  after  deviating  from  the  mission  trajectory  (Lilley,  2011).  Major 
Adams  was  killed  in  the  resulting  breakup,  and  the  X-15  program  lost  its  funding  the  next 
fiscal  year.  NASA  reports  the  lessons  from  the  X-15  were  the  basis  of  the  flight  control 
system  for  the  Apollo  Lunar  Excursion  Module,  and  heavily  influenced  the  Space  Shuttle 
program  (Lilley,  2011). 

Concurrently  to  the  X-15  project,  adaptive  controls  were  being  proposed  and  put 
to  the  test  in  many  other  fields,  such  as  spacecraft  control  and  machine  tooling.  In 
partnership  with  the  Bendix  Corporation,  the  USAF  experimented  with  utilizing  adaptive 
controls  with  feedback  for  metal  cutting,  one  of  the  first  such  attempts  (Ulsoy  &  Koren, 
1989).  The  goal  was  to  account  for  tool  wear  in  the  machining  process,  both  by 
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measurement  and  by  position  feedback.  At  the  same  time  NASA  was  investigating  the 
possibility  of  developing  adaptive  controls  for  space  vehicles,  though  the  stated  purpose 
was  to  extend  the  lifetime  of  spacecraft  by  minimizing  the  fuel  consumption  required  to 
perfonn  the  mission  (Taylor,  1969).  Gain  scheduling  was  heavily  favored  in  these 
applications,  as  the  required  computing  power  to  perform  nonlinear  adaptive  control  was 
yet  to  come.  Gain  scheduling  requires  significant  up  front  processing  power,  in  order  to 
develop  all  the  linear  controls  needed,  which  also  suffered  from  the  relative  lack  of 
computing  power  of  the  era. 

These  projects  led  naturally  to  the  inclusion  of  adaptive  controls  in  the 
development  of  the  International  Space  Station,  to  assist  with  stability  during  the 
construction  phase  and  overall  attitude  control  (Ih,  Wang,  &  Leondes,  1985).  In  this  we 
see  the  emergence  of  developing  inherently  robust  controls  across  the  entire  mission 
profile,  where  only  the  parameters  change  between  the  different  phases.  Also  emerging 
at  this  point  are  many  different  attempts  to  limit  uncertainty  for  space  based  systems, 
such  as  the  use  of  Kalman  filters,  model  reference  controls,  and  linear  quadratics  (Govin, 
Claudinon,  &  Lanninat,  1981).  Indeed,  the  realization  that  slight  disturbances  can  lead  to 
system  failure  has  forced  the  development  these  inherently  robust  controls,  and  leads 
developers  to  calculate  uncertainty  and  feed  it  into  the  control  parameters  adjustment 
during  operations  (Kosut,  1989).  As  late  as  1983,  the  equations  required  to  perform 
adaptive  control  for  complex  space  structures  were  too  complicated  to  be  run  in  real  time 
(Schaechter,  1983). 

B.  RECENT  DEVELOPMENT 

Development  up  to  the  1990s  left  space  operations  needing  a  controller  that 
functions  well  across  a  system  with  changing  operational  conditions.  Serious  study  at 
this  point  focused  on  utilizing  the  dominant  control  method  of  industry,  the  proportional 
integral  derivative  (PID)  controller,  which  was  developed  for  naval  autopilots  (Xu,  Shum, 
Lee,  &  Kanade,  1991).  The  adaptive  PID  controller  works  for  a  wide  array  of  processes 
without  requiring  significant  characterization  of  the  system,  all  while  outperforming  other 
control  methods,  such  as  the  generalized  predictive  control  system  (Ho,  Rad,  Chan,  & 
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Wong,  1999).  This  is  due  in  part  to  the  flexibility  of  the  controller,  which  is 
demonstrated  in  the  adaptive  control  scheme  used  by  this  simulation.  Due  to  the  nature 
of  the  PID  equation,  developers  can  take  specific  portions  of  the  controller  and  adapt 
them  for  a  specific  use,  such  as  the  Proportional  Derivative  controllers  proposeded  for 
use  in  space  manipulators  (Ehrenwald  &  Guelman,  1998).  PID  and  PID  derivative 
adaptive  controls  show  significant  promise  when  applied  to  space  based  controls. 
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III.  SIMULATION  METHODOLOGY 


A.  INTRODUCTION 

The  purpose  of  the  simulation  is  to  answer  our  research  question  “Are  accurate 
and  stable  DTC  possible  for  satellites  in  orbit?”  Simulating  ground  conditions  allows  for 
continuation  and  validation  of  any  results  in  the  laboratory.  The  satellite  was  modeled  in 
SIMULINK  and  the  satellite  block  diagram  is  as  shown  in  Figure  1 . 


Figure  1.  Spacecraft  Attitude  Control  Simulation. 

The  system  receives  user  commanded  Euler  angles  and  sends  them  into  the 
trajectory  generator.  Trajectory  is  then  sent  from  the  generator  to  the  Actuators  and 
Control  block.  Inside  of  the  control  portion  is  the  ability  for  the  system  to  select  between 
the  PID  controller  and  the  modified  PID  controller,  and  to  enable  or  disable  feed  forward 
control.  Leaving  the  control  block  is  the  product  of  the  control  selections  and  the 
feedback  controls,  the  commanded  torques,  which  are  run  through  the  CMG  actuators, 
Dynamics  block,  and  then  Kinematics  block.  The  resulting  spacecraft  Euler  angles  are 
finally  output  from  the  simulation.  Disturbances  due  to  gravity  are  then  fed  back  into  the 
Dynamics  block.  Feedback  is  fed  through  the  Sensors  and  Observers  block,  where  the 
system  can  select  between  utilizing  ideal,  sensor,  or  observer  feedback  and  enable  the 
Kalman  and/or  low  pass  filters.  The  simulation  utilizes  initial  conditions  akin  to  those 
expected  in  a  laboratory  via  the  model  initialization  function,  which  can  be  found  in 
Appendix  A. 
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B. 


MODIFIED  PID  CONTROLLER 


The  concept  of  the  PID  controller  comes  from  the  Astrom  and  Wittenmark 
definition: 


1  r  /  \  „  de  ^ 


!(0  =  KC  e(t)  +  y\e(s)ds  +  T1 


dt 


(1) 


V  n  o  "V 

where  u(7)  is  the  process  or  control  input,  e{t)  is  the  error  defined  as  e  =  uc  -y  where 


uc  is  the  reference  value,  y{t )  is  the  process  output,  and  Kc,  2} ,  and  TD  are  constants 

used  to  tune  the  controller  (Astrom  &  Wittenmark,  2008,  p.  376).  The  controller  sends 
torque  commands  seeking  the  desired  angle  to  achieve  the  commanded  trajectory.  As  the 
torque  command  is  and  acceleration,  the  controller  seeks  the  correct  acceleration  to 
produce  the  desired  angle.  From  Equation  1  we  see  that  the  proportional  portion 
depends  on  current  error,  the  integral  portion  accumulates  past  error,  and  the  derivative 
portion  extrapolates  future  error.  The  modified  PID  controller  is  adapted  to  accept  the 
selected  type  of  feedback  from  the  feedback  block,  and  outputs  the  control  signal.  By 
passing  the  derivative  action  channel  of  the  control  to  a  separate  channel,  the  modified 
PID  controller  avoids  differentiating  noisy  signals,  increasing  control  efficiency.  The 
modified  PID  block  diagram  is  as  shown  in  Figure  2. 


Figure  2.  Modified  PID  Controller. 
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c. 


STEERING 


The  spacecraft  achieves  attitude  control  via  CMG  steering,  within  the  CMG 
Steering  subsystem,  which  is  contained  in  the  Actuators  and  Control  subsystem  as  shown 
in  Figure  1.  The  simulation  uses  a  minimum,  non-redundant  CMG  array  of  three  CMGs 
and  a  balance  mass  to  offset  gravity  gradient  disturbances.  The  CMG  Steering  subsystem 
utilizes  the  Moore-Penrose  Pseudoinverse  Steering  Law  as  described  by  Bong  Wie  (Wie, 
2008,  pp.  465-466).  Wie’s  description  assumes  a  CMG  array  as  shown  in  Figure  3. 


The  CMG  array  produces  torque  by  manipulating  the  angular  momentum  of  the  CMGs 
by  changing  their  angular  velocity.  The  angular  momentum  vector  is  defined  as: 


-  cos  (3  sin  8X 

-  cos  8 2 

cos  P  sin  83 

cos  84 

cosSl 

+ 

-  cos  P  sin  8 2 

+ 

-  COS  C>3 

+ 

cos  P  sin  8a 

sin  j3  sin  8X 

sin  p  sin  S2 

sin  p  sin  S3 

sin  p  sin  84 
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where  P  is  the  pyramid  skew  angle,  {//.}  is  the  angular  momentum  vector  of  the  i  th 
CMG,  and  {4}  is  the  gimbal  angle  of  the  i  th  CMG.  Taking  the  time  derivative  of  the 
angular  momentum  vector  results  in  the  equation: 

{*}-£{*} 


=Zl«,(!t)Hl  (J) 

1=1 

where  |c>}  is  the  gimbal  angle  vector,  {a,.}  is  the  i th  column  of  [^4] ,  and  [ A ]  is  the 

Jacobian  matrix  defined  as: 

—  cos  [3  cos  8y  sin  S2  cos  [3  cos  S3  -sin£4 

[A]-  -sin£j  -  cos  [3  cos  S2  sin  £3  cos/Jcosfii,  .  (4) 

sin  [3  cos  8X  sin  /?  cos  S2  sin  p  cos  c>3  sin  p  cos  S4 

For  a  commanded  control  torque  input,  {w},  the  CMG  momentum  rate  command  j h  j  is 

defined  as: 

[h}  =  -{u}-{co}x{h\  (5) 

where  ico }  is  the  angular  velocity.  Finally  the  Moore-Penrose  Pseudoinverse  Steering 
Law,  also  known  as  pseudoinverse  steering  logic,  may  be  obtained  as: 


where  \A]+  is  the  pseudoinverse  of  [ A ]  if: 


ifi- infix 


A+=Ar(AAry 


AA+A  =  A 


A+AA+  =  fi4 
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K) 

=  AA+ 

(10) 

'  =A+ A 

(11) 

where  A  =  [A\  .  Generally  speaking,  the  pseudoinverse  [A]  ,  is  the  conjugate  transpose 
of  [^4] ,  a  relationship  otherwise  known  as  Hermitian  matrices.  As  the  simulation  uses 
only  3  CMGs,  [A]  in  the  simulation  is  a  3x3  square,  real  matrix.  Therefore,  the 
conjugate  transpose  of  [A]  is  equal  to  the  inverse  of  [d]  and  the  pseudoinverse  of  [ A ] : 

A*  =  A~1  =  A+  (12) 

for  the  simulation’s  specific  case  only.  The  subsystem  takes  the  controller  commanded 
torque  and  removes  the  current  torque,  giving  the  actual  torque  command.  This  is 

multiplied  with  [A]+  to  find  the  commanded  gimbal  rate,  which  is  in  turn  sent  to  the 
gimbal  actuators.  Integrating  the  gimbal  rate  gives  the  gimbal  state,  which  is  an  input  of 
[  A ] .  The  subsystem  outputs  angular  momentum  rate  or  torque  rate  after  multiplying  [d] 
with  the  gimbal  rate.  The  CMG  Steering  subsystem  is  as  shown  in  Figure  4. 


Figure  4.  CMG  Steering. 


D.  DYNAMICS 

The  Dynamics  section  was  derived  from  the  law  of  conservation  of  angular 
momentum,  and  the  definition  of  angular  momentum: 

{/f,}  =  [/]{®}  (13) 
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where  j  Hs }  is  the  spacecraft  angular  momentum,  [/]  is  the  moment  of  inertia  matrix, 
and  {&>}  the  spacecraft  angular  velocity.  From  the  definition  of  angular  momentum 
comes  Euler’s  second  law  of  motion: 


{• », } 


d{H,} 


dt 


(14) 


where  { /W /v }  is  the  sum  of  the  external  moments  of  force,  or  torques  around  an  axis  in 
the  inertial  frame.  Euler’s  second  law  as  applied  to  rotating  frames  is: 


{MRor}  =  d^!r-+{M} 


dt 


(15) 


where  {Mrot}  is  the  sum  of  the  torques  in  the  rotational  frame  and  { M  ] ,  the  rotation  of 
the  body  frame,  is  defined  as  (Reeves,  1999,  pp.  324-325): 

{M\  =  {a\x{H,}.  (16) 

From  change  in  angular  momentum  we  use  Equation  13  to  find  change  in  angular 
velocity  and  angular  velocity,  outputting  both.  The  spacecraft  dynamics  block  diagram  is 
as  shown  in  Figure  5. 


Imo*  omegadot= 


Hs  = 


N  B 


Figure  5.  Spacecraft  Dynamics. 
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E.  KINEMATICS 

The  Kinematics  section  begins  with  the  angular  velocity  of  the  spacecraft, 
(cox,coy,coz) ,  and  adds  the  angular  velocity  of  the  spacecraft’s  orbit,  resulting  in 

(cdx,co2,co2)  ,  the  angular  velocity  with  respect  to  orbit.  This  angular  velocity  is 

transformed  into  a  quaternion,  which  is  a  specific  type  of  complex  number,  by  the 
relationship: 

[e]=i{n}®[e]  d7) 

where  [Q\  is  the  quaternion,  is  the  rate  of  change  of  the  quaternion,  and  {Q}  is  the 

angular  velocity  with  respect  to  orbit.  Expanded,  Equation  17  becomes  the  kinematic 
differential  equation  for  quaternions: 

qx  0  -CO*,  cdx  qx 

q2  1  -co 3  0  cox  a>2  q2 

q2  2  (Dr,  -cox  0  co2  q2 

qA  -cox  -co2  -co2  0  g4 

where  qx  is  the  scalar  part,  the  values  q2 ,  q3 ,  q4  make  up  a  vector  part,  and  [  ]  is  a 

matrix  multiplication  skew-symmetric  version  of  the  cross  product  in  Equation  17  (Wie, 
2008,  p.  426).  One  benefit  of  using  Equation  18  to  find  the  quaternion  is  that  a 
quaternion  calculated  via  this  method  will  always  be  nonnalized  (Wie,  2008,  p.423). 
Once  the  quaternion  \Q\  is  calculated,  the  direction  cosine  matrix  (DCM)  is  calculated 
via  the  relationships: 

M'=[eim[e]‘  a?) 

[ef  =  ~ q 1  (20) 

-q2 

_“^4_ 
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where  {P}  is  the  point  subject  to  rotation  as  described  by  quaternion  [Q],  '  is  the 


resulting  point,  and  \Q\  is  the  conjugate  quaternion  (MathWorks,  2013).  When  Equation 


19  is  expanded,  the  following  DCM  can  be  extracted  (Wie,  2008,  p.  335): 

1-2  (q]+ql)  2(qtq2+q3q4)  2  (qxq3~q2q4) 

[DCM]=  2(yq2q,-q3q4)  1-2  (q[ +q")  2(q2q3+qiq4) 

2(q3q\+q2q4)  2{q3q2-qxq4 )  l-2(q*+ql) 


(21) 


The  quaternion  block  is  required  in  order  to  avoid  singularities  during  the  simulation  that 
result  when  using  Euler  angles  alone.  The  Apollo  program’s  inertial  measurement  unit 
was  highly  susceptible  to  gimbal  lock  thanks  to  using  Euler  angles  without  quaternions 
and  three  axis  gimbals,  leading  to  dangerous  losses  of  position  accuracy  during  several 
missions  (Hoag,  1963).  The  block  diagram  of  the  quaternion  and  DCM  calculations  is  as 
shown  in  Figure  6. 
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Figure  6. 


Quaternion  and  DCM  Calculation. 
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The  DCM  is  also  defined  by  writing  out  the  result  of  the  three  body  axis  rotations: 

c(0)c(y/)  c(0)s  (y/)  -s  (9) 

[ DCM]=  s(^)s(6i)c(^)-c(^)s(^)  s(^)s(6*)s(^)-c(^)c(^)  s(^)c(6*)  (22) 
_c(^)s(0)c(^)  +  s(^)s(^)  c(^)s(6>)s(^)-s(^)c(^)  c(/)c(0)_ 

where  c(x)  =  cos(x) ,  ,s(x)  =  sin  (x) ,  and  assuming  the  coordinate  transfonnation 

[Ct]WHC,}(e)^[Cr](y)  where  [C,](x)  is  the  rotation  matrix  about  the  i  th  axis 

with  an  angle  of  x  (Wie,  2008,  p.  328).  Having  used  Equation  21  to  populate  the  DCM 
with  values,  the  simulation  uses  Equation  22  to  find  the  corresponding  Euler  angles. 
Pitch,  or  9 ,  can  be  found  by  taking  the  value  of  the  first  row,  third  column  of  the  DCM, 
DCMU ,  applying  the  corresponding  portion  of  Equation  22,  and  then  isolating  6 : 

DCMu=-sin(9) 

-sin  1  (DCM13)  =  -shT1  (-sin(0)) 

9  = -sin  l(DCMn).  (23) 

With  9  found,  the  remainder  of  the  Euler  angles  can  be  similarly  isolated  and  solved: 


<f>  =  sin 


DCM23  ' 
cos (9)  y 


(24) 


i//  =  sin 


DCM[2  ' 
cos(6*) 


(25) 


where  DCM  is  the  element  in  row  x ,  column  y  of  the  DCM.  The  block  diagram  of 
the  Euler  angle  generator  is  as  shown  in  Figure  7. 
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C12 


Figure  7.  Euler  Angle  Generator. 

After  calculating  the  Euler  angles,  they  are  output  from  the  Simulation.  The 
block  diagram  of  the  Kinematics  section  is  as  shown  in  Figure  8. 


Spacecraft  Kinematics 


N  O 


Euler 


*CD 


Euler(deg) 


Figure  8.  Spacecraft  Kinematics. 


F.  DISTURBANCES 

The  only  disturbances  modeled  were  the  gravity  gradient  torque,  or  the  torque  due 
to  unequal  gravitational  pull  on  the  disparate  parts  of  the  satellite.  This  is  the  only 
disturbance  expected  to  be  observed  in  a  laboratory  follow  up  of  this  simulation.  Gravity 
gradient  torque  is  defined  as: 

P6) 

Ko 

where  j T „  j  is  the  gravity  gradient  torque,  //is  the  Earth’s  gravitational  coefficient,  R0 
is  the  distance  from  Earth’s  center,  { u\  is  the  unit  vector  towards  Nadir,  and  [/]  is  the 
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spacecraft  moment  of  inertia  matrix  (Reeves,  1999,  p.  324).  The  unit  vector  towards 
Nadir  is  picked  out  of  the  DCM: 


DCMn 

DCM2i 

DCM}3 


(27) 


The  gravity  gradient  torque  is  output  to  the  Dynamics  Block,  and  the  block  diagram  is  as 
shown  in  Figure  9. 


Gravity  Gradient  Torques 


Figure  9.  Gravity  Gradient  Torques. 

G.  SENSORS  AND  OBSERVERS 

The  enabling  concept  behind  the  feedback  capabilities  of  the  simulation  is 
contained  in  the  Sensors  and  Observers  block.  Euler  angles  and  angular  rates  enter, 
combined  in  a  state  vector,  which  can  be  shunted  along  different  paths,  from  an  unedited 
ideal  full  state  feedback  to  sensor  feedback,  observers,  state  feedback,  Kalman  filtered 
sensors,  and  Lowpass  filtered  sensors.  Once  the  simulation  selects  anything  other  than 
ideal  full  state  feedback,  a  sensor  is  simulated.  The  sensor  is  the  ideal  feedback  signal 
with  a  random  number  generator  adding  noise.  Next  in  the  process  flow  is  the  Observer 
Subsystem,  containing  a  Luenberger  observer  and  a  PID  observer.  The  purpose  of  any 
observer  is  to  create  a  complete  state  vector  from  incomplete  and/or  noisy  observations. 
The  PID  observer  relies  upon  the  same  concept  as  the  PID  controller,  or  measuring 
current  error,  accumulating  past  error,  and  predicting  future  error.  The  PID  observer 
equation,  derived  from  Equation  1,  is: 
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(28) 


j#j  =  |  {(£>}  dt 

=  JJ  ( Kpe  +  dt  +  KDe^  dtdt 

where  j#j  is  the  estimated  Euler  angle,  {cb}  is  the  estimated  angular  velocity,  e(t  j  is  the 

error  defined  as  e  =  uc-y  where  uc  is  the  reference  value,  y  (/)  is  the  process  output, 
and  Kp  ,  K j ,  and  are  gain  constants  used  to  tune  the  observer.  The  PID  controller 
operates  by  accepting  angle  data  which  is  used  to  estimate  the  angular  acceleration  used 
to  produce  the  angle  data.  Integrating  this  estimate  gives  estimated  angular  velocity  {&>}  , 

integrating  again  produces  estimated  Euler  angles  j<9j .  Tuning  the  gain  constants  Kp , 
K ! ,  and  KD  such  that  the  estimates  approach  the  actual  values,  or  {&>}—>{&>}  and 
|(9|  — »  {#}  ,  achieves  what  is  known  as  full  state  feedback. 

The  other  choice  for  observer,  the  Luenberger  observer,  also  has  the  ability  to 
receive  the  control  signal  from  the  Dynamics  block,  labeled  TotalControl  on  the  block 
diagram.  The  Luenberger  observer  operates  on  a  system  defined  as: 


x(t)  =  Ax(t}  +  Bu{t} 

(29) 

y(t)  =  Cx(t) 

(30) 

where  x ( / )  is  the  process  state,  u\t)  is  the  process  inputs,  y(t  )  is  the  process  outputs, 

and  A,B ,  C,  and  D  are  gain  constants  for  tuning  the  observer  (Luenberger,  1979). 
Assuming  such  a  system,  the  Luenberger  observer  is  defined  as: 

z(t)  =  y4z(7)  +  [f](  y(7)-Cz(7))  +  fh/(f)  (31) 

e(t)  =z{t)-x(t) 

=  (A-[L\C)(z{t)-x{t ))  (32) 

=  (A-[L\C)e(t) 

Where  z{t)  is  the  estimated  process  state,  e(?)is  the  observer  error,  and  [f]  is  the 
observer  gain  matrix.  The  observer  gain  matrix  is  tuned  such  that  the  observer  error 
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converges  to  zero.  The  simulated  Luenberger  observer  takes  in  the  Euler  angles,  finds 
the  observer  error  compared  to  the  process  state,  and  uses  Equation  32  to  output  full  state 
feedback.  The  key  difference  between  the  PID  and  Luenberger  observer  is  the 
Luenberger’s  ability  to  sum  the  derivative  portion  of  the  estimation  with  the  velocity 
component.  By  doing  this,  the  Luenberger  does  not  put  the  derivative  portion  through  the 
additional  differentiation  used  by  the  PID  observer,  as  shown  in  Equation  28.  By 
manipulating  the  derivative  portion  less,  the  Luenberger  observer  propagates  less  noise 
than  the  PID  observer.  As  implemented  in  the  simulation,  the  Luenberger  observer 
equation  is: 

{<?}  =  {  ({  (Kpe  +  K,\  (e)  dt)  +  KDe) .  (34) 

The  simulated  Luenberger  observer  is  enhanced  to  accept  feedback  just  like  the  PID 
controller  and  features  the  TotalControl  option  as  previously  discussed.  The  block 
diagram  of  the  PID  and  Luenberger  observers  are  as  shown  in  Ligure  10. 


Ligure  10.  Observer  Subsystem. 


Immediately  after  the  Observer  block  is  the  ability  to  enable  a  Kalman  filter. 
Kalman  filters  are  a  means  of  isolating  a  system  state  buried  inside  an  otherwise  noisy 
signal.  The  Kalman  filter  operates  on  a  state  defined  as: 

W = [5]{**-i} + Mk-J + k-J  (35) 


21 


where  { xk }  is  the  state,  [5]  is  the  transition  matrix,  [l>]  is  the  control  matrix,  { uk }  is 
the  control,  and  { wk }  is  the  process  noise,  all  at  time  k  (Welch  &  Bishop,  2006).  The 
filter  makes  observations  defined  as: 

{z*MFlW+k}  (36) 

where  { zk }  is  the  observation,  [F]  is  the  observation  matrix,  and  { vk }  is  the  observation 

noise.  Both  process  and  observation  noise  are  assumed  to  normally  distributed  with 
covariance  Q  and  R  respectively.  With  these  definitions,  the  Kalman  filter  makes 
estimates  via  this  equation: 

{xt)  =  Kk{zk}  +  (  l-X,){xH}  (37) 

where  { xk }  is  the  estimated  state  and  Kk  is  the  Kalman  gain.  Using  these  relationships  a 

Kalman  filter  makes  a  prediction  of  the  state  and  covariance,  calculates  gain,  and  updates 
its  state  and  covariance  model. 

In  the  simulation,  the  state  vector  is  broken  apart,  and  the  Euler  angles  are  run 
separately  from  the  angular  velocity  through  corresponding  Kalman  filters.  This  enables 
the  simulation  to  find  the  true  error  in  the  Kalman  values  by  comparing  the  filtered  Euler 
angles  to  be  compared  to  the  true  Euler  angles.  The  Kalman  Filter  block  diagram  is  as 
shown  in  Figure  11. 


Figure  1 1 .  Kalman  Filter. 
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The  final  piece  of  the  Sensors  and  Observers  block  is  the  Lowpass  filter,  which 
attempts  to  eliminate  extraneous  high  frequency  noise  in  our  signal: 

;l=Jy^+2^M+i  (38) 

s  K  +  2C/M+1 

where  s  is  the  signal,  5 '  is  the  filtered  signal,  and  Q)p ,  C p ,  0)_ ,  C,z  are  parameters  tuned 

to  the  desired  filter  level  (Wie,  2008,  p.  137).  When  the  lowpass  filter  is  selected,  the 
simulation  breaks  out  each  of  the  six  components  of  the  state  vector,  utilizes  individual, 
separately  tunable  filters,  and  pushes  the  filtered  signal  to  the  Controller  block.  The 
Lowpass  filter  block  diagram  is  as  shown  in  Figure  12. 


CL> 


1  /wz*  2s?+2*dampZ/wzsM 

1/wp*2s?+2*dampP/wpsM 

1  /wz*  2s?+2*dampZ/wzsH 

1  /wp*  2s?+2*damp  P/wpsM 

1  /wz*  2s?+2*dampZ/wzsM 

1  /wpA  2s?+2*damp  P/wpsM 

1  /wz1'  2s?+2*dampZ/wzsM 

1  /wpA  2s?+2’damp  P/wpsH 

1  /wz*  2  s£ +2  *d  a  mp  Z/wzs+ 1 

1/wpA2s?+2*dampP/wps+1 

1  /wz*  2s?+2*dampZ/wzs+1 

1  /wp*  2s?+2*damp  P/wpsH 

KjD 

u  (filtered) 


Figure  12.  Lowpass  Filter. 

The  Sensors  and  Observers  block  offers  a  wide  array  of  abilities  for  the 
simulation.  After  passing  through  the  enabled  blocks,  the  full  state  feedback  is  sent  to  the 
Dynamics  block.  The  Sensors  and  Observers  block  diagram  is  as  shown  in  Figure  13. 
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Figure  13.  Spacecraft  Sensors  and  Observers. 


H.  CHAPTER  SUMMARY 

This  chapter  covered  the  simulation  as  it  was  written,  including  the  many 
different  possible  configurations.  Each  section  included  the  founding  principle  and 
associated  equations  that  led  the  development  of  the  respective  block. 
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IV.  ANALYSIS  AND  RESULTS 


A.  INTRODUCTION 

To  answer  the  research  question  completely  requires  that  the  DTC  are  effective, 
accurate,  and  stable.  This  requires  at  least  two  different  methods  for  analysis.  The 
selected  methods  are  the  Monte  Carlo  analysis  and  phase  portrait  analysis.  Both  methods 
will  demonstrate  that  DTC  are  effective,  or  that  the  simulated  DTC  do  in  fact  function. 
The  Monte  Carlo  analysis  will  also  demonstrate  the  effects  on  accuracy  due  to  DTC.  The 
phase  portrait  analysis  will  demonstrate  the  stability  of  DTC. 

B.  MONTE  CARLO  ANALYSIS 

The  Monte  Carlo  method,  created  by  Stanislaw  Ulam,  draws  its  name  from 
Ulam’s  uncle,  who  frequented  the  Monte  Carlo  casino  (Metropolis,  1987).  The  premise 
is  that  by  running  simulations  numerous  times  and  recording  the  results,  the  outcome  can 
be  accurately  characterized.  This  is  akin  to  visiting  a  casino  many  times  to  detennine  the 
odds  of  any  game  of  chance. 

In  order  to  profile  the  outcome  of  the  DTC  in  the  simulation,  a  Matlab  m-file  was 
written,  iteratively  calling  the  simulation.  Each  iteration  was  distinguished  by  increasing 
the  mass  loss  of  the  spacecraft  by  0.1  percent,  beginning  at  zero  percent  and  finishing  at 
90  percent  mass  loss.  Mass  loss  was  simulated  by  reducing  the  simulated  spacecraft 
moment  of  inertia  matrix.  This  assumes  that  in  the  case  of  sudden  mass  loss,  the  satellite 
remains  functional.  Each  iteration,  the  spacecraft  was  commanded  to  perfonn  a  30 
degree  yaw,  and  the  tracking  error  was  recorded.  In  addition  to  the  mean  error,  the 
standard  deviation  of  the  error  was  also  recorded.  The  m  file  was  written  as  shown  in 
Appendix  B. 

The  Monte  Carlo  method  applied  to  tracking  error  does  an  excellent  job  of 
characterizing  the  accuracy  of  the  DTC  when  a  spacecraft  mass  is  suddenly  lost,  and  by 
returning  these  results,  shows  that  the  DTC  of  the  simulation  do  function  correctly.  The 
DTC  never  experienced  greater  than  a  mean  0.9  degree  tracking  error,  with  a  standard 
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deviation  that  was  always  less  than  1.1  degrees.  The  results  are  as  shown  in  Figure  14, 


where 


is  the  mean  measured  error,  and 


is  the  standard  deviation. 


Figure  14.  Monte  Carlo  Analysis  Results. 


C.  PHASE  PORTRAIT  ANALYSIS 

The  premise  behind  phase  portrait  analysis  is  to  characterize  the  behavior  of  a 
system  with  regard  to  stability  of  equilibrium  points.  In  this  case,  the  behavior  of  the 
spacecraft  DTC  was  examined  around  the  commanded  state  vector.  In  order  to 
demonstrate  stability,  the  axes  of  angular  rate  and  angular  position  in  the  inertial  frame 
were  chosen.  If  the  motion  described  circles  into  the  commanded  point,  then  stability  has 
been  achieved. 

To  create  phase  portraits,  a  Matlab  m-file  was  written  which  iteratively  calls  the 

simulation.  Each  portrait  reflects  a  fixed  mass  loss,  with  many  different  values  of  initial 
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position  and  initial  velocity  given.  Mass  loss  was  simulated  by  reducing  the  inertia 
matrix.  This  assumes  that  in  the  case  of  sudden  mass  loss,  the  satellite  remains 
functional.  Each  iteration  changes  the  initial  position  and  initial  velocity,  then  records 
the  results  as  the  simulation  is  run  for  a  set  time  of  100s.  The  omega  value  was  already 
calculated  by  the  Dynamics  section,  as  variable  w.  In  order  to  record  the  desired  results, 
a  new  variable  had  to  be  captured,  the  position  in  the  inertial  frame.  This  was 
accomplished  by  adding  an  integrator  block  to  the  angular  speed  in  the  Dynamics  section, 
and  capturing  the  angle  as  the  variable  wAngles,  as  shown  in  Figure  15. 


M=H  X  w 


Figure  15.  Omega  Angle  Production. 

For  each  iteration  the  spacecraft  was  commanded  to  perfonn  a  30  degree  yaw, 
with  the  omega  angles  and  omega  value  recorded,  then  plotted  against  each  other  in 
graphs  for  each  of  the  inertial  frame  dimensions.  The  m-file  was  written  as  shown  in 
Appendix  C. 

The  generated  phase  trajectory  portraits  successfully  characterizes  the  stability  of 
the  DTC  when  spacecraft  mass  is  suddenly  lost  and  the  spacecraft  is  imparted  with 
sudden  changes  in  state.  Stability  is  demonstrated  by  the  boundedness  of  the  phase 
trajectories  (Slotine  &  Fi,  1991).  Stable  trajectories  are  generally  expected  to  approach 
the  origin,  unstable  trajectories  generally  diverge  from  the  origin.  Examples  by  Slotine 
and  Fi  of  phase  trajectories  with  differing  stability  can  be  found  in  Figure  16. 
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curve  1  -  asymptotically  stable 
curve  2  -  marginally  stable 
curve  3  -  unstable 


Figure  16.  Phase  Portrait  Stability. 


By  returning  these  results,  the  portraits  demonstrates  that  the  DTC  of  the  simulation  do 
function  correctly.  Each  set  of  phase  portraits  is  grouped  by  the  inertial  frame  dimension 
it  comes  from,  and  the  portraits  show  either  zero,  45,  or  90  percent  mass  loss  of  the 
spacecraft.  The  phase  portraits  are  shown  in  Figures  17  through  25. 
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Figure  17.  0%  Mass  Loss  X  Phase  Portrait. 
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Figure  18.  45%  Mass  Loss  X  Phase  Portrait. 
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Figure  19.  90%  Mass  Loss  X  Phase  Portrait. 
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Figure  20.  0%  Mass  Loss  Y  Phase  Portrait. 


Angle(deg) 


Figure  2 1 .  45%  Mass  Loss  Y  Phase  Portrait. 
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Figure  22.  90%  Mass  Loss  Y  Phase  Portrait. 
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Figure  23.  0%  Mass  Loss  Z  Phase  Portrait. 
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Figure  24.  45%  Mass  Loss  Z  Phase  Portrait. 
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Figure  25.  90%  Mass  Loss  Z  Phase  Portrait. 
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D.  CHAPTER  SUMMARY 

This  chapter  covered  the  analysis  methods  used  to  answer  the  research  question. 
The  Monte  Carlo  method  proved  effective  in  characterizing  the  tracking  error  of  the 
system  when  subjected  to  mass  loss.  The  phase  portraits  provide  insight  into  the  stability 
of  the  system  when  subjected  to  sudden  mass  loss,  change  in  velocity,  and  change  in 
position.  In  both  examples  the  fewest  possible  changes  were  made  to  the  base  simulation 
to  protect  the  integrity  of  the  results.  The  changes  were  limited  to  commenting  out,  or 
removing,  initial  conditions  that  interfered  with  data  collection,  and  adding  a  single 
integrator  and  variable  output  set  for  phase  portrait  generation. 
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V.  CONCLUSIONS 


A.  INTRODUCTION 

The  development  of  adaptive  controls  and  the  enabling  technology  for  them  has 
reached  a  point  where  new  and  innovative  uses  are  now  becoming  possible.  By 
introducing  a  modified  adaptive  PID  controller  with  adaptive  feed  forward  control  to  this 
simulated  spacecraft,  it  is  demonstrated  that  the  controls  have  achieved  significant 
damage  tolerance.  This  comes  at  a  time  where  the  need  for  such  tools  has  never  been 
greater,  and  the  dangers  associated  with  space  operations  grow  daily. 

B.  ADAPTIVE  CONTROL 

Space  operations  are  going  to  become  more  difficult  and  dangerous  in  the 
foreseeable  future.  Significant  action  is  required  to  protect  the  space  assets  that  the 
United  States  has  become  dependent  on.  Adaptive  control  holds  the  key  for  many  of  the 
solutions  currently  proposed  to  mitigate  these  risks. 

The  damage  tolerant,  adaptive  controller  as  implemented  in  the  simulation 
provides  adaptive  control  with  accuracy  and  stability.  The  accuracy  is  evidenced  by  the 
minimal  tracking  error  of  the  system  even  when  approaching  a  90  percent  mass  loss.  The 
stability  is  evidenced  by  the  phase  portraits  in  each  inertial  frame  recovering  from  the 
sudden  change  in  state  even  when  the  change  was  accompanied  by  a  loss  of  mass.  The 
research  question  was:  Are  DTC  possible  for  satellites  in  orbit?  Yes,  the  satellite 
simulation  achieved  accurate  and  stable  DTC. 

C.  APPLICATION  OP  STUDY 

Introducing  DTC  like  those  simulated  here  is  achievable  with  today’s  technology. 
As  the  largest  limiting  factor  historically  has  been  processing  power,  the  feasibility  of 
retrofitting  operational  spacecraft  with  DTC  is  limited,  but  not  necessarily  impossible. 
Modern  spacecraft  should  have  no  technical  limitations  prohibiting  the  implementation  of 
DTC.  Writing  DTC  into  spacecraft  controllers  may  even  lower  the  cost  and  risk  of  the 
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program  as  a  whole  (Wertz  &  Larson,  1999).  DTC  clearly  exhibit  a  robust  ability  for  all 
manner  of  spacecraft,  warranting  further  study  and  operational  adoption. 


D.  AREAS  TO  CONDUCT  FURTHER  RESEARCH 

The  logical  next  step  is  to  validate  the  simulation  results  in  a  laboratory.  The 
simulation  was  designed  to  enable  laboratory  testing,  by  including  the  Gravity  Gradient 
disturbances  alone  and  setting  the  altitude  to  zero,  as  would  be  observed  in  the  laboratory. 
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APPENDIX  A.  MODEL  INITIALIZATION  FUNCTION 


%  Spacecraft  Dynamics  Simulation 
%clear  all;close  all;clc; 
warning  off  all; 

%  Simulation  run  parameters 
Rate=34.3;  DeltaT=l/Rate; 

%  Maneuver  Parameters 

roll=0;  pitch=0;  yaw=30;  roll=roll*pi/180;  pitch=pitch*pi/180; 
yaw=yaw*pi / 180; 

SlewTime=5;  SACengage=0; 

%ConstantsO 

Re=6371 . 2e3;mu=398601 .  2e9; 
constant 

%Spacecraft  orbit 
h=0  ; 

R=Re+h; 

we=0. 000072921158553; 
sec (Vallado) 
wo=sqrt (mu/ (Re+h) A3) ; 
incln=36.6001*pi/180; 
epsilon=12*pi/180;  alphao=0; 
uo=0;  nuo=0;  %Start  S/C  beneath  subsolar  point 

beta=60;  gamma=1.5; 

a=0 . 545491852;  b=0 . 314939867;  c=0 . 704226952;  %Assumed  spacecraft 
rectangular  size 

Area=[b*c  a*c  a*b] ;  %projected  area~mA2  in  body  x,y,z 

directions 

density=4 . 39e-14 ; 

kpre=-9 . 9639/24/3600/180*pi*0;  %nodal  precession  constant  assumed 

zero  here 

wn=kpre* (Re/ (Re+h) ) A3 . 5*cos (incln) ;  %nodal  precession  (zero 
eccentricity) 

V=wo* (Re+h) ; 

rho=asin (Re/ (h+Re) ) ;  %earth  angular  radius 

Cd=2.5;  psun=4.5E-6;  %Drag  coefficient  and  solar  pressure 
constant~N/mA2 

Kaero=-0 . 5*Cd*VA2 ;  Psolar=2*psun;  %constants  for  aero  and  solar  torque 
calculation 

dL= [0.002  0.002  0.008];  %predicted  distance  between  cp  and  eg 

Kme=2 ,3390e-005; 

%Spacecraft  Magnetic  Properties  (assumed) 

mresid=[0  0  0.01];  %Spacecraft  residual  magnetic  moment 

M=mresid;  %Magnetic  unit  dipole  vector 

K=7 . 943el5; 


%earth  radius  and  universal  gravitation 


%orbit  altitude 

%orbit  radius  from  center  of  earth 
%earth's  angular  velocity  rad/solar 

%orbit  angular  velocity 

%Inclination  assumed  Monterey  Latitude 
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%Spacecraft  Inertia  conditions  (ACTUAL  values ....  not  the  ones  assumed 
in  the  feedforward  control  calculation) 

Ix=90;  Iy=100;  Iz=250  ; 

Ixy=-10;  Iyz=20;  Ixz=-10; 

Imo=[Ix  -Ixy  -Ixz; 

-Ixy  Iy  -Iyz ; 

-Ixz  -Iyz  Iz]  ;  %Moment  of  inertia  matrix 

Iinv=inv ( Imo) ;  %Moment  of  inertia  inverse  goes  in  dynamics 

block 

%Spacecraft  initial  Euler  state  angles  and  rates 
phio=0 ; thetao=0 ; psio=0 ;  %Initial  Euler  Angles 
phidoto=0 ; thetadoto=0 ; psidoto=0 ;  %Initial  Euler  Rates 

%Calculation  of  initial  quaternion  (qo)  and  angular  momentum  (Ho) 

sl=sin(phio/2) ;s2=sin (thetao/2 ) ;s3=sin (psio/2) ;cl=cos (phio/2) ;c2=cos (th 

etao/2) ;c3=cos (psio/2) ; 

qlo=sl*c2*c3-cl*s2*s3; 

q2o=cl*s2*c3+sl*c2*s3;  %Wie  pg .  321 

q3o=cl*c2*s3-sl*s2*c3; 

q4o=cl*c2*c3+sl*s2*s3; 

Sl=sin(phio) ;S2=sin (thetao) ;S3=sin (psio) ; Cl=cos (phio) ; C2=cos (thetao) ; C3 
=cos (psio) ; 

wxo=phidoto-psidoto*S2-wo*S3*C2 ; 

wyo=thetadoto*Cl+psidoto*C2*Sl-wo* (C3*C1+S3*S2*S1 ) ; 
wzo=psidoto*C2*Cl-thetadoto*Sl-wo* (S3*S2*C1-C3*S1) ; 
qo=[qlo  q2o  q3o  q4o] ; 

Ho=Imo* [wxo  wyo  wzo]  '  ; 
norm (Ho) *  1000 ; 

%Calculate  eclipse  time  for  comparison  with  EPS  calculations 
Te=100 . 87*2*V/2/pi; 

%CMG  Properties  (in  degrees) 

%beta=90*pi/180;  %Skew  angle  in  degrees  converted  to 

radians 

beta= [ 90 ; 90 ; -90 ; 0 ] ;  beta=beta . *pi . / 180 ; 

Gimbal0= [-30*pi/180;  90*pi/180;  -30*pi/180; 0] ;  %  Initial  Gimbal  angles 
for  0  H  spin  up 

w  wheel=2800* (2*pi/60 ) ;  %Wheel  speed  in  RPM  converted 

to  rad/s 

lwheel=0 . 0614*1 . 3558179483314;  %  Wheel  inertia  in  slug-ftA2 

converted  (exact)  to  kg  mA2 

h  wheel=Iwheel*w  wheel;  %  CMG  Wheel  Angular  Momentum 

%  [Fossen] 's  adaptive  feedforward  parameters 
ETA=-100 ;  LAMBDA=0 . 5 ;  uffGain=l; 

%FEEDFORWARD  control  based  on  "assumed"  Spacecraft  Inertia  conditions 
I x= 119. 1259;  Iy=150.6615;  Iz=106.0288;  %  ACTUAL  VALUES  ARE  Ix=90; 
Iy=100;  Iz=250  ; 

Ixy=-15 . 7678;  Iyz=22 . 31637;  Ixz=-6 . 54859;  %  ACTUAL  VALUES  ARE  Ixy=-10; 
Iyz=20;  Ixz=-10; 

THETAo= [ lx  Ixy  Ixz  Iy  Iyz  Iz  0  0  0]; 


36 


%  FEEDBACK  CONTROL 

%PDI  Controller  Gains  %  TUNED  Well  for  Presence  of  LOWPASS  (Kp=0.5; 
Kd=Kp*750;  Ki=0.1) 

%Kp=0.5;  Kd=Kp*750;  Ki=0.1; 

%PDI  Controller  Gains  %  TUNED  WELL  FOR  NO-NOISE  (Kp=l;  Kd=2000*Kp; 
Ki=5 ) 

Kp  •••  1  ;  Kd=Kp*6000;  Ki=Kp*75; 

%PDI  Controller  Gains  %  TUNED  WELL  FOR  NOISE  (Kp=l;  Kd=Kp*3000;  Ki=l) 
%Kp=l ;  %Kd=Kp*3000 ;  %Ki=l; 

%PID  Controller  Gains  tuned  well  for  uff  and  ufb  decoupling 
Kpx=20;  Kdx=500;  Kix=0.1; 

%  Noise  Parameters 

NoiseVariance=le- 9;  BADNoiseVariance=NoiseVariance*le3 ; 

%Bandpass  filter  pole  per  Bong  Wie  2nd  Edition  pg  137-138 
wp=10*pi/SlewTime;  %  Product  of  pole  frequency  and  zero  frequency 
establishes  max  phase  lag 
wz=2*wp; 

dampZ=l;  %  >0  but  small  for  good  tracking.. 
dampP=dampZ ; 

%  Luenberger  Observer  Gains 

Maxl=300;  lambdal=12.5  ;  lambda2=50  ;  lambda3=  200; 

Kpo=MaxI* (lambdal* ( lambda2+ lambda 3 ) +Iambda2*lambda3) /10; 

Kdo=MaxI* ( lambdal +lambda2+l ambda3 ) /10; 

Kio=MaxI* lambdal *lambda2*lambda3/ 9; 

%  PID  Observer  Gains 
KpoPID=30000 ; 

KdoPID=3000 ; 

KioPID=30 ; 
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APPENDIX  B.  DAMAGE  TOLERANT  CONTROLS  MONTE  CARLO 

CODE 


%Iterative  Mass  Loss  Analysis  Scott  Nakatani 

%This  file  calls  the  Spacecraft  simulation  with  increasing  mass  loss  to 
%assess  the  ability  of  the  adaptive  controls  to  compensate  for  battle 
%damage,  characterized  as  a  percent  mass  lost 
clear  all;  close  all;  clc 

StepSize=0 . 001;  %ENTER  STEP  SIZE  HERE  where  1  =  100% 

Steps=round (0 . 9/StepSize) ;  %Calculates  #  of  steps  to  take 

Loss=[];  PercentLoss=0 ;  MeanErrorT= [ ]  ;  StdDevT=[];  DeltaT=0 . 001  ; 


%Actual  Spacecraft  Inertia  conditions: 

%These  are  separate  from  the  simulation's  feedforward  assumed  values 
Ix=90;  Iy=100;  Iz=250  ;  Ixy=-10;  Iyz=20;  Ixz=-10; 

ImoI=[Ix  -Ixy  -Ixz;  -Ixy  Iy  -Iyz;  -Ixz  -Iyz  Iz] ;  %Moment  of 
inertia 

for  ii  =  1: Steps 

%Sends  %  complete  status  to  command  window 
PercentCompleted=100* (ii/Steps) 

%INVERSE  of  Moment  of  inertia  matrix 
Imo= (1-PercentLoss) . *ImoI; 

%Accumulates  values  of  battle  damage  percent 
Loss=[Loss;  PercentLoss ] ; 

%CALLS  THE  SIMULATION 

sim(  ' IterativeSpacecraf tModel . mdl '  ); 

%Accumulate  Tracking  Error  &  Standard  Deviation 
MeanErrorT= [MeanErrorT;  norm (MeanError, 2 ) ] ; 

StdDevT= [ StdDevT;  norm ( StdDev, 2 ) ] ; 

%Iterates  mass  lost 
PercentLoss=StepSize*ii ; 

end 


%Insure  array  dimensions  match 
Loss=Loss (1 :max (size (StdDevT) ) ) ; 

MeanErrorT=Loss (1 :max (size (StdDevT) ) )  ; 

%Plot  the  results  of  Monte  Carlo  analysis 
figure ( 1 ) ; 

plot (Loss* 100 , MeanErrorT, ' linewidth 3 ) ;  grid  on;  ylabel (' Tracking 
Error  (Degrees) fontsize 16) ;xlabel ( 'Mass  Loss 
(Percent) ' , ' fontsize ' , 16) ; 

hold  on;  plot (Loss*100 , StdDevT, ; linewidth ', 3 ) ; 

legend ('ll \muAo I  I  ' ,  '  I  I \sigmaAo II',' fontsize ' , 16,  ' location '  ,  'Northwest ' ) 

r 

hold  off 
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APPENDIX  C.  PHASE  PORTRAIT  CODE 


%Stability  Phase  Portrait  Generator  Scott  Nakatani 

%This  file  calls  the  Spacecraft  simulation  with  increasing  angular 
change 

%to  assess  the  stability  of  the  adaptive  control  solution  via  creating 
a 

%phase  portrait 

clear  all;  close  all;  clc 

Omega= [ ] ;  OmegaAngles= [ ] ; 

OmegaX= [ ] ;  OmegaY= [ ] ;  OmegaZ= [ ] ; 

OmegaAnglesX= [ ] ;  OmegaAnglesY= [ ] ;  OmegaAnglesZ= [ ] ; 

%Actual  Spacecraft  Inertia  conditions: 

%These  are  separate  from  the  simulation's  feedforward  assumed  values 
Ix=90;  Iy=100;  Iz=250  ;  Ixy=-10;  Iyz=20;  Ixz=-10; 

Imo=[Ix  -Ixy  -Ixz;  -Ixy  Iy  -Iyz;  -Ixz  -Iyz  Iz] ;  %Moment  of  inertia 

Imo=Imo . * ( . 1) ;  %INPUT  MOMENTUM  LOSS  HERE  where  1  =  100%  mass,  0%  loss 
Steps=4;  %Steps  are  in  all  quadrants,  total  steps  =  2*StepsA2 
StepSize=2*pi/360;  %Stepsize  is  input  in  radians 

for  ii  =  1: Steps 

%Iterates  initial  positive  omega  values  &  subordinate  conditions 

omegaOX= (ii) * (StepSize)  ; 

omegaOY= (ii) * (StepSize) ; 

omegaOZ= (ii) * (StepSize) ; 

omega0= [omegaOX,  omegaOY,  omegaOZ]'; 

H0=lmo*omega0 ; 

for  iii  =  l:Steps 

%Sends  %  complete  status  to  command  window 
PercentCompleted= (50* ( ii- 1 ) /Steps ) + (50/ Steps* iii/ Steps ) 

%Iterates  initial  positive  omega  angles 

wAngleX= (iii) * (StepSize) ; 

wAngleY= (iii) * (StepSize) ; 

wAngleZ= (iii) * (StepSize)  ; 

wAngle0= [wAngleX,  wAngleY,  wAngleZ]; 

%CALLS  THE  SIMULATION 

sim(  ' PhasePortraitSpacecraf tModel . mdl '  ); 

%Accumulate  resulting  omega  values  and  angles 
OmegaX= [OmegaX  w ( : , 1) ] ; 

Omega Y= [ Omega Y  w ( : , 2 ) ] ; 

OmegaZ= [OmegaZ  w(:,3)]; 

OmegaAnglesX= [OmegaAnglesX  wAngles ( : , 1) ] ; 

OmegaAnglesY= [OmegaAnglesY  wAngles ( : , 2 ) ] ; 

OmegaAnglesZ= [OmegaAnglesZ  wAngles ( : , 3) ] ; 
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%Iterates  initial  negative  omega  angles 

wAngleX= ( — iii ) * (StepSize) ; 

wAngleY= (-iii) * (StepSize) ; 

wAngleZ= (-iii) * (StepSize)  ; 

wAngleO= [wAngleX,  wAngleY,  wAngleZ]; 

%CALLS  THE  SIMULATION 

sim(  ' PhasePortraitSpacecraf tModel . mdl '  ); 

%Accumulate  resulting  omega  values  and  angles 
OmegaX= [OmegaX  w ( : , 1)  ]  ; 

Omega Y= [ Omega Y  w  (  :  ,  2 )  ]  ; 

OmegaZ= [OmegaZ  w(:,3)]; 

OmegaAnglesX= [OmegaAnglesX  wAngles (:  ,1)  ]  ; 

OmegaAnglesY= [OmegaAnglesY  wAngles ( : , 2 ) ] ; 

OmegaAnglesZ= [OmegaAnglesZ  wAngles (:  ,3)  ]  ; 

end 

end 

for  iiii  =  1: Steps 

%Iterates  initial  negative  omega  values  &  subordinate  conditions 

omegaOX= (-iiii) * (StepSize)  ; 

omegaOY= (-iiii) * (StepSize)  ; 

omegaOZ= (-iiii) * (StepSize)  ; 

omegaO= [omegaOX,  omegaOY,  omegaOZ]'; 

HO=Imo*omegaO ; 

for  iiiii  =  1: Steps 

%Sends  %  complete  status  to  command  window 

PercentCompleted=50+ (50* (iiii-1) /Steps) t (50/Steps*iiiii/Steps) 

%Iterates  initial  positive  omega  angles 
wAngleX= (iiiii) * (StepSize) ; 
wAngleY= (iiiii) * (StepSize)  ; 
wAngleZ= (iiiii) * (StepSize)  ; 
wAngleO= [wAngleX,  wAngleY,  wAngleZ]; 

%CALLS  THE  SIMULATION 

sim(  ' PhasePortraitSpacecraf tModel . mdl '  ); 

%Accumulate  resulting  omega  values  and  angles 
OmegaX= [ OmegaX  w ( : , 1 ) ] ; 

OmegaY= [ Omega Y  w ( : , 2) ] ; 

OmegaZ= [OmegaZ  w(:,3)]; 

OmegaAnglesX= [OmegaAnglesX  wAngles (: ,1) ] ; 

OmegaAnglesY= [OmegaAnglesY  wAngles ( : , 2 ) ] ; 

OmegaAnglesZ= [OmegaAnglesZ  wAngles ( : , 3) ] ; 

%Iterates  initial  negative  omega  angles 
wAngleX= (-iiiii) * (StepSize) ; 
wAngleY= (-iiiii) * (StepSize)  ; 
wAngleZ= (-iiiii) * (StepSize) ; 
wAngleO= [wAngleX,  wAngleY,  wAngleZ]; 
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%CALLS  THE  SIMULATION 

sim(  ' PhasePortraitSpacecraf tModel . mdl '  ); 


%Accumulate  resulting  omega  values  and  angles 
OmegaX= [ OmegaX  w ( : , 1 ) ] ; 

Omega Y= [ Omega Y  w ( : , 2 ) ] ; 

OmegaZ= [OmegaZ  w(:,3)]; 

OmegaAnglesX= [OmegaAnglesX  wAngles (  :  ,1)  ]  ; 
OmegaAnglesY= [OmegaAnglesY  wAngles ( : , 2)  ]  ; 
OmegaAnglesZ= [OmegaAnglesZ  wAngles ( :  ,  3)  ]  ; 

end 

end 


%Converts  omega  angles  &  values  to  degrees 
OmegaX=OmegaX . *360/ (2*pi) ; 

OmegaY=OmegaY . *360/ (2*pi) ; 

OmegaZ=OmegaZ . *360/ (2*pi)  ; 
OmegaAnglesX=OmegaAnglesX . *360/ (2*pi)  ; 
OmegaAnglesY=OmegaAnglesY . *360/ (2*pi)  ; 
OmegaAnglesZ=OmegaAnglesZ . *360/ (2*pi)  ; 


%Plots  the  phase  portraits 
figure ( 1 )  ; 

plot (OmegaAnglesX, OmegaX, ' linewidth ' , 1 ) ;  grid; 
xlabel ( ' Angle_x (deg) ' , ' fontsize '  ,  16)  ;  ; 
ylabel ( ' \omega_x (deg/s) ' , ' fontsize ' , 16) ; 
figure  (2 ) ; 

plot (OmegaAnglesY, OmegaY, ' linewidth ', 1 ) ;  grid; 
xlabel ( ' Angle_y (deg) ' , ' fontsize ' , 16) ; ; 
ylabel ( ' \omega_y (deg/s) ' , ' fontsize ' , 16) ; 
figure  (3); 

plot (OmegaAnglesZ , OmegaZ , ' linewidth ', 1 ) ;  grid; 
xlabel ( 'Angle_z (deg) ' , ' fontsize ' , 16) ; ; 
ylabel ( ' \omega_z (deg/s) ' , ' fontsize ' , 16) ; 
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